library(raster)
## Warning: package 'raster' was built under R version 3.5.3
## Loading required package: sp
library(sp)
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.5.3
## rgdal: version: 1.4-2, (SVN revision 814)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/ASUS PC/Documents/R/win-library/3.5/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/ASUS PC/Documents/R/win-library/3.5/rgdal/proj
## Linking to sp version: 1.3-1
library(wesanderson)
## Warning: package 'wesanderson' was built under R version 3.5.3
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 3.5.2
library(tmap)
## Warning: package 'tmap' was built under R version 3.5.2
b1 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF1.tif')
b2 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF2.tif')
b3 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF3.tif')
b4 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF4.tif')
b5 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF5.tif')
b6 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF6.tif')
b7 <- raster('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF7.tif')
par(mfrow = c(2,2))
plot(b2, main = "Blue", col = gray(0:100/100))
plot(b3, main = "Green", col = gray(0:100/100))
plot(b4, main = "Red", col = gray(0:100/100))
plot(b5, main = "NIR", col = gray(0:100/100))
stack false colour and rgb
landsatRGB <- stack(b3, b2, b1)
landsatFCC <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
plotRGB(landsatFCC, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
## full landsat
files <- paste0('C:/Users/ASUS PC/Documents/R/r_remote/1981/TIF', 1:7, ".tif")
landsat_1981 <- stack(files)
names(landsat_1981) <- c('blue','green','red','nir','swir1','tir','swir2')
extent(landsat_1981)
## class : Extent
## xmin : 252769.2
## xmax : 288859.2
## ymin : 522889.4
## ymax : 558259.4
landsat_81_sub <- subset(landsat_1981, 1:3)
nlayers(landsat_81_sub)
## [1] 3
pairs(landsat_1981[[1:2]],main='blue v. green')
pairs(landsat_1981[[3:4]],main='red v. nir')
### normalised index
nd <- function(img, a, b){
ba <- img[[a]]
bb <- img[[b]]
index <- (ba-bb)/(ba+bb)
return(index)
}
### bare soil index. a= swir1, b= red, c= nir, d= blue
bli <- function(img,a,b,c,d){
ba <- img[[a]]
bb <- img[[b]]
bc <- img[[c]]
bd <- img[[d]]
index <- ((ba+b3)-(bc+bd))/((ba+b3)+(bc+bd))
}
ndvi <- nd(landsat_1981, 4,3)
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
ndwi <- nd(landsat_1981, 2,4)
plot(ndwi, col = rev(terrain.colors(10)), main = 'Landsat-NDWI')
ndbi <- nd(landsat_1981, 5,4)
plot(ndbi, col = rev(terrain.colors(10)), main = 'Landsat-NDBI')
ndbl <- bli(landsat_1981, 5,3,4,1)
plot(ndbl, col = rev(terrain.colors(10)), main = 'Landsat-NDBLI')
hist(ndvi,
main = "Distribution of NDVI values",
xlab = "NDVI",
ylab= "Frequency",
col = "green",
xlim = c(-1, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-1,1, 0.05), labels = seq(-1,1, 0.05))
hist(ndwi,
main = "Distribution of NDWI values",
xlab = "NDWI",
ylab= "Frequency",
col = 'steelblue3',
xlim = c(-1, 0.8),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-1,0.8, 0.05), labels = seq(-1,0.8, 0.05))
hist(ndbi,
main = "Distribution of NDBI values",
xlab = "NDBI",
ylab= "Frequency",
col = "navajowhite",
xlim = c(-1, 0.7),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-1,0.7, 0.05), labels = seq(-1,0.7, 0.05))
summary(ndvi)
## layer
## Min. -4.857143e-01
## 1st Qu. 4.214876e-01
## Median 5.370370e-01
## 3rd Qu. 6.056338e-01
## Max. 7.815126e-01
## NA's 8.130520e+05
summary(ndwi)
## layer
## Min. -6.666667e-01
## 1st Qu. -4.918033e-01
## Median -4.368932e-01
## 3rd Qu. -3.543307e-01
## Max. 6.279070e-01
## NA's 8.130520e+05
summary(ndbi)
## layer
## Min. -9.736842e-01
## 1st Qu. -2.456140e-01
## Median -1.834320e-01
## 3rd Qu. -5.633803e-02
## Max. 5.902778e-01
## NA's 8.134540e+05
veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main = 'Veg cover')
water <- reclassify(ndwi, cbind(-Inf, -0.3, NA))
plot(water, main = 'Water cover')
built <- reclassify(ndbi, cbind(-Inf, -0.2, NA))
plot(built, main = 'Built cover')
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Veg Over RGB")
plot(veg, add=TRUE, legend=FALSE)
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Water Over RGB")
plot(water, add=TRUE, legend=FALSE)
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Built Over RGB")
plot(built, add=TRUE, legend=FALSE)
writeRaster(veg, filename= "1981_veg-mask.tif", overwrite=TRUE)
writeRaster(water, filename= "1981_water-mask.tif", overwrite=TRUE)
writeRaster(built, filename = "1981_built-mask.tif", overwrite=TRUE)